Модели временных рядов с гетероскедастичностью
ARCH модель
Пусть нам дан временной ряд, удовлетворяющий модели авторегрессии порядка \(p\) \[x_{t}=c+\phi_1x_{t-1}+...+\phi_{p}x_{t-p}+u_{t}\] где \(u_{t}\)-белый шум \[Eu_{t} = 0\] \[Eu_{t}u_{s} = σ^2\] при \(s=t\) иначе 0
Оптимальным линейным прогнозом такого ряда является условное математическое ожидание \(x_{t}\) при условии, что нам известна история процесса до момента \(t-1\) включительно \[E[x_{t}|x_{t-1},x_{t-2,....}]=c+\phi_1x_{t-1}+...+\phi_{p}x_{t-p}\] (1)
Если предположить стационарность процесса \(x_{t}\) то безусловное математическое ожидание процесса не зависит от времени и равно \[E[x_{t}] = c+\phi_1E[x_{t-1}]+...+\phi_{p}E[x_{t-p}]= c+\phi_1E[x_{t}]+...+\phi_{p}E[x_{t}]\] Отсюда
\[E[x_{t}]=(с/(1-\phi_1-...-\phi_{p}))\]
Однако из выражения (1) видно, что условное математическое ожидание от времени зависит.
На практике довольно часто встречается необходимость прогнозирования не только средних значений, но и дисперсии. Предположим, что и сам процесс шума \(u_{t}\)(точнее его квадрат) в свою очередь следует процессу авторегрессии поряда \(m\)
\[u_{t}^2=\eta+\alpha_1u_{t-1}^2+...+\alpha_m u_{t-m}^2+v_{t}\] (2) Где \(v_{t}\)- процесс белого (но уже другого) шума \[Ev_{t} = 0\] \[Ev_{t}v_{s} = \lambda^2\] при \(s=t\) иначе 0
При прогнозировании ряда \(x_{t}\) мы показали, что \(u_{t}\)- ошибка прогнозирования. Из представления (2) следует, что оптимальным линейным прогнозом \(u_{t}^2\) будет в свою очередь следующее условное математическое ожидание.\[E[u^2_{t}|u_{t-1}^2,u^2_{t-2},....]=\eta+\alpha_1u^2_{t-1}+...+\alpha_{m}u^2_{t-m}\] (3) Определение. Процесс белого шума, удовлетворяющий (2) называют авторегрессионным гетероскедастическим процессом порядка m. Обозначают это следующим образом \(u_{t}\sim ARCH(m)\). Такой класс процессов впервые был предложен Engle в 1982 году. В 2003 году вклад Engle был отмечен Нобелевской премией по экономике. Так как процесс \(u_{t}\) - случайный, a квадрат его \(u_{t}^2\) естественно должен быть всегда неотрицательным, то от представления (2) естественно потребовать, что при любых реализациях \(u_{t}\) - его квадрат был всегда нетрицательным. Этого можно добится следеющими предположениями 1. Процесс \(v_{t}\) ограничен снизу \(-\eta\), где \(\eta>0\) 2. Все \(\alpha_{i}>0;i=1,...,m\) Естественно также предполагать стационарность процесса \(u_{t}^2\). Это достигается предположением, что корни соответсвующего характеристического уравнения многочлена \[1-\alpha_1z-...-\alpha_{m}z^{m}=0\] должны лежать вне единичного круга \(|z|≤1\)
Это предположение, а вместе с предположением 2. эквивалентны тому, что \[\alpha_1+...+\alpha_{m}<1\]
При выполнении всех выше сделанных предположений безусловную дисперсию процесса \(u_{t}\) можно вычислить по следующей формуле \[Eu_{t}^2=(\eta/(1-(\alpha_1+...+\alpha_{m})))\] Существует альтернативное и более удобное представление ARCH(m). Предположим, что
\[u_{t}=e_{t}\sqrt(h_{t})\], (4)
где \(e_{t}\) -последовательность независимых одинаково распределенных случайных величин с нулевым математическим ожиданием и единичной дисперсией
\[E[e_{t}] = 0\], \[D[e_{t}] = 1\],
а \(h(t)\) удовлетворяют соотношению \[h_{t}=\eta+\alpha_1u^2_{t-1}+...+\alpha_{m}u^2_{t-m}\] (5) Из соотношения (4) следует, что \[E[u_{t}^2|u_{t-1},u_{t-2},...]=\eta+\alpha_1u^2_{t-1}+...+\alpha_{m}u^2_{t-m}\]
следовательно \(u_{t}\sim ARCH(m)\)
Оценивание параметров ARCH модели.
Предположим, что у нас есть линейная регрессионная модель с шумом, который в свою очередь предствляется в виде ARCH процеса.\[y_{t}=x_{t}′\beta+u_{t}\]
где \(x_{t }\)-вектор из независимых (предикторов) переменных, среди которых могут быть и так называемые лаговые переменные от \(y_{t}\), например \(y_{t-1},y_{t-2}\) и т.д. А \(u_{t}\) удовлетворяет соотношениям (4) и (5). Удобно занумеровать значения имеющихся ряда наблюдений \(y_{t}\) индексами, начиная с \(t=-m+1,-m+2,...,0,1, ...,T\). Первые \(m\) значенией мы будем использовать для построения условного распредения, а значения с индексами от 1 до \(Т\) непосредстевнно для оценивания. Таким образом для произвольного момента времени \(t\) обозначим через \(Y_{t}\) следующий вектор \[Y_{t}=(y_{t},y_{t-1},...,y_0,...,y_{-m+1},x_{t}′,x_{t-1},...,x′_0,...,x_{-m+1}′\]
Пусть \(e_{t}\)- последовательность независимых стандартно нормально распределенных случайных величин независимых одновременно от \(x_{t}\) и от \(Y_{t-1}\). Тогда, при сделанных выше предположениях условное распределение \(y_{t}\) при условии, что \(x_{t}\) и \(Y_{t-1}\) известно будет нормальным с математическим ожиданием \(x_{t}′\beta\) и дисперсией \(h_{t}\) \[f(y_{t}|x_{t},Y_{t-1})=(1/(\sqrt{(2\pi h_{t}})))exp(-(((y_{t}-x_{t}′\beta)^2)/(2h_{t})))\] где \[h_{t}=\mu+\alpha_1(y_{t-1}-x_{t-1}′\beta)^2+...+\alpha_{m}(y_{t-m}-x_{t-m}′\beta)^2=\] (8) Запишем (8) в виде cкалярного произведения двух векторов \[=[z_{t}(\beta)]′\delta\] где \[\delta=(\mu,\alpha_1,...,\alpha_{m})′\],
а \[[z_{t}(\beta)]′=[1,(y_{t-1}-x_{t-1}′\beta)^2,...,(y_{t-m}-x_{t-m}′\beta)^2]\]
Обозначим через \(\theta\) вектор параметров подлежащих оценке. Это будет вектор вида
\[\theta=[\beta′,\delta′]\]
Тогда логарифм условной относительно первых m наблюдений функции правдоподобия запишется следующим образом \[L(\theta) = \sum_{t=1}^{T}ln(f(y_{t}|x_{t},Y_{t-1};\theta))= -(T/2)ln(2\pi)-(1/2)\sum_{t=1}^{T}ln h_{t}-(1/2)\sum_{t=1}^{T}(((t_{t}-x_{t}′\beta)^2)/(h_{t}))\] Для получения оценок максимального правдоподобия чаше всего используют какой-нибудь градиентный метод. Познакомимся с ARCH(m) процессом в следующим фрейме. Здесь моделируется ARCH(2) с добавленным ARMA(1,1) процессом.
Моделируемая модель выглядит так: \[y_t=\mu + \phi_1y_{t-1}+\theta_1a_{t-1}+a_t\], где \(a_t=\sqrt{h_t}\epsilon_t\)$ \[h_t= \omega+\alpha_1h_{t-1}+\alpha_2h_{t-2}\]
R код моделирования выглядит так. Для работы ARCH процессами нам понадобится библиотеки rugarch и fBasics
library(rugarch)## Loading required package: parallel library(fBasics)## Loading required package: timeDate## Loading required package: timeSeries## ## Rmetrics Package fBasics## Analysing Markets and calculating Basic Statistics## Copyright (C) 2005-2014 Rmetrics Association Zurich## Educational Software for Financial Engineering and Computational Science## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.## https://www.rmetrics.org --- Mail to: info@rmetrics.org r = rnorm(1000)
spec <- ugarchspec(variance.model = list(model ="sGARCH",garchOrder = c(2,0)), mean.model = list(armaOrder = c(1, 1)),distribution.model = "norm")
fit <- ugarchfit(r, spec = spec)
mu <- 0.5
AR <- 0.4
MA <- 0.7
omega <- 3
alpha1 <- 0.45
alpha2 <- 0.3
fit@fit$ipars[1]<-mu #mu
fit@fit$ipars[2]<-AR #ar1
fit@fit$ipars[3]<-MA #ma1
fit@fit$ipars[4]<-0.0 #arfima
fit@fit$ipars[5]<-0.0 #arch
fit@fit$ipars[7]<-omega #omega
fit@fit$ipars[8]<-alpha1 #alpha1
fit@fit$ipars[9]<-alpha2 #beta1
#fit@fit$coef <- fit@fit$ipars
# fit@fit$matcoef <-fit@fit$ipars
sim <- ugarchsim(fit,n.sim=500, n.start=1, m.sim=1, startMethod="sample")
retval <- sim@simulation$seriesSim
matplot(retval,type = "l",lwd = 3,col = "blue",main = "Simulated time series")
Вот так выглядит гистограмма смоделированного ряда по сравнению с плотностью нормального распределения
histPlot(as.timeSeries(sim@simulation$seriesSim))
Cмоделированный процесс для дисперсии \(\sigma_t^2\) можно посмотреть вот так:
sigma_t <- sim@simulation$sigmaSim
matplot(sigma_t,type = "l",lwd = 3,col = "blue",main = "Simulated sigma processs")
Результаты оценивания методом максимального правдоподобия
fitarch <- ugarchfit(retval, spec = spec)
fitarch@fit$coef## mu ar1 ma1 omega alpha1 alpha2
## 0.8818839 0.3977756 0.7124691 2.8144422 0.4654690 0.2833069Остатки после удаления модели \(\epsilon_t\).
resid <- fitarch@fit$residuals
matplot(resid,type = "l",lwd = 3,col = "blue",main = "Residuals")
QQ-график остатков
qqnorm(resid ,lwd=3,col = "blue")
GARCH модель
В предыдущем разделе \(h(t)\) удовлетворяло соотношению \[h_{t}=\mu+\alpha_1u^2_{t-1}+...+\alpha_{m}u^2_{t-m}\]
Более общим будет процесс в котором условная дисперсия зависит не от конечного, а бесконечного числа задержек \(u^2_{t-j}, j=1,2,...\) В операторном виде это можно записать так. \[h_{t} = \mu+\pi(B)u_{t}^2 (1)\] \[\pi(B) = \sum_{j=1}^{\infty}\pi_{j}B^{j}\]
Эта идея реализуется, когда \(\pi(B)\) представим в виде отношения двух конечных полиномов \[\pi(B)=((\alpha(B))/(1-\delta(B)))=((\alpha_1B+...+\alpha_{m}B^{m})/(1-\delta_1B-...-\delta_{r}B^{r}))\] здесь мы предположим, что корни \(1-\delta(B)\) лежат вне единичного круга \(|B|≤1\)
Умножив обе части (1) на \(1-\delta(B)\) получим \[[1-\delta(B)]h_{t}=[1-\delta(B)]\mu+\alpha(B)u_{t}^2\] или \[h_{t} = κ+\delta_1h_{t-1}+...+\delta_{r}h_{t-r}+\alpha_1u_{t-1}^2+...+\alpha_{m}u_{t-m}^2 (2)\] где \(κ = [1-\delta_1-...-\delta_{r}]\mu\) Соотношения (2) и называют процессом обобщенной условной авторегрессии с гетероскедастичностью. (GARCH(r,m) - модель). Модель впервые была предложена Bollerslev в 1986 году.
Неотрицательность для \(u_{t}^2\) обеспечивается предположением, что \(κ≥0,\alpha_{i}≥0,\delta_{i}≥0; i=1,...,p=max(r,m)\). Здесь мы доопределили нулями недостающие либо \(\alpha_{i}\) либо \(\delta_{i}\). Cтационарность гарантируется тем, корни характеристического моногочлена вида \[1+(\alpha_1+\delta_1)z+...+(\alpha_{p}+\delta_{p})z^{p}=0 \] лежат вне единичного круга \(|z|≤1\)
Совместно с предположением о нетрицательности \(\alpha_{i}\) и \(\delta_{i}\) предположение о стационарности эквивалентно тому, что \[(\alpha_1+\delta_1)+...+(\alpha_{p}+\delta_{p})<1 \]
При сделанных предположениях безусловное математическое ожидание для \(u_{t}^2\) будет равно \[Eu_{t}^2=(κ/(1-((\alpha_1+\delta_1)+...+(\alpha_{p}+\delta_{p}))))\]
Ниже приведен пример кода на R для моделирования GARCH(2,2) модели с добавленной ARMA(1,1). Обратите внимание мы поменяем распределение шума \(\epsilon_t\) на стандартное распределение Стьюдента с нулевым математическим ожиданием и дисперсией 1.
r = rnorm(1000)
spec <- ugarchspec(variance.model = list(model ="sGARCH",garchOrder = c(2,2)), mean.model = list(armaOrder = c(1, 1)),distribution.model = "std")
fit <- ugarchfit(r, spec = spec)## Warning in arima(data, order = c(modelinc[2], 0, modelinc[3]), include.mean
## = modelinc[1], : possible convergence problem: optim gave code = 1 mu <- 0.5
AR <- 0.4
MA <- 0.7
omega <- 3
alpha1 <- 0.25
alpha2 <- 0.3
betta1 <- 0.2
betta2 <- 0.2
fit@fit$ipars[1]<-mu #mu
fit@fit$ipars[2]<-AR #ar1
fit@fit$ipars[3]<-MA #ma1
fit@fit$ipars[4]<-0.0 #arfima
fit@fit$ipars[5]<-0.0 #arch
fit@fit$ipars[7]<-omega #omega
fit@fit$ipars[8]<-alpha1 #alpha1
fit@fit$ipars[9]<-alpha2 #beta1
fit@fit$ipars[10]<-betta1 #betta1
fit@fit$ipars[11]<-betta2 #betta2
sim <- ugarchsim(fit,n.sim=500, n.start=1, m.sim=1, startMethod="sample")
retval <- sim@simulation$seriesSim
matplot(retval,type = "l",lwd = 3,col = "blue",main = "Simulated time series")
Обратите вниманиe как изменилась гистограмма смоделированного ряда по сравнению с плотностью нормального распределения
histPlot(as.timeSeries(sim@simulation$seriesSim))
Оценку неизвестных параметров GARCH модели проводится аналогично как и для ARCH методом максимального правдоподобия. Соответсвующее выражение для логарифма функции правдоподобия можно взять из работы Bollerslev 1986 года.
fitarch <- ugarchfit(retval, spec = spec)
fitarch@fit$coef## mu ar1 ma1 omega alpha1
## 0.551175043 0.349609974 0.699684796 2.150324985 0.296499259
## alpha2 beta1 beta2 shape
## 0.159025512 0.000905739 0.500935799 99.996990477Подробнее ознакомитсяя с GARCH(m,k) процессом можно в следующим фрейме. Здесь моделируется GARCH(1,1) с добавленным ARMA(1,1) процессом и нормальным распределением шума \(\epsilon_t\).
Прогнозирование GARCH моделей
Прогнозирование рассмотрим на примере модели AR(1)$GARCH(1,1).
\[x_t=\phi_1x_{t-1}+\theta_1\sigma_{t-1}\epsilon_{t-1}+\sigma_t\epsilon_t\] \[\sigma_t^2= \omega+\alpha a_t^2+ \beta \sigma^2_{t-1}\] \[a_t=\sigma_t\epsilon_t; \epsilon_t\sim N(0,1)\]
Обозначим для момента времени \(t\) прогноз на \(l\) моментов времени вперед \(\widetilde{x}_t(l)\). Как и в ARMA моделях условное математическое ожидание будет минимальным в среднеквадратичном прогнозом
\[\widetilde x_t(l) = E(x_{t+l}|x_1,x_2,...,x_t)=\phi_1\widetilde x_{t}(l-1)\], где \(\widetilde x_t(k) = E(x_{t+k}|x_1,x_2,...,x_t)\) if \(k>0\)
и
\(\widetilde x_t(k) = x_{t+k}\) для \(k\le 0\).
Пусть \(\widetilde \sigma_t^2(l)\) прогноз квадрата волатильности в момент времени \(t\) на \(l\) моментов вперед. \[\widetilde \sigma_t^2(l) = E(\sigma^2_{t+l}|\sigma^2_1,\sigma^2,...,\sigma^2_t)=\alpha\widetilde a^2_{t}(l-1)+\beta\widetilde \sigma^2_{t}(l-1)\], где
\(\widetilde a^2_t(k) = 0\) if \(k>0\);
иначе
\(\widetilde a^2_t(k) = x_{t+k}-\widetilde x_{t}(k)\) if \(k \le 0\).
Обратим внимание, что прогноз самого ряда \(x_t\) точно такой же, как если бы у нас просто была модель AR(1). GARCH отражается только на доверительных интервалах прогноза Следующий фрейм дает представление о возможностях прогнозирования по GARCH моделям
Во многом по этим причинам позднее появились обобщения модели для учета GARCH характера модели в прогнозах самого значения процесса \(x_t\), а не только в доверительных интервалов. Исторически первой появилась модель
GARCH-M модель
В финансовых рядах доходность актива может зависить от волатильности. Рассмотрим модель GARCH-M символ “М” в названии модели означает “in the maean”. Простейшая GARCH-M(1,1) модель записывается следующим образом
\[r_t=\mu+c\sigma_t^2+a_t,a = \sigma_t\epsilon_t\] \[\sigma_t^2= \alpha_0+\alpha_1a^2_{t-1}+\beta_1\sigma_{t-1}^2\]
Знак параметра \(c\) определяет тип зависимости динамики рейта от волатильности (положительную и отрицательную корреляцию между этими параметрами). Озакомиться с процессом ,увидеть влияние параметров на процесс можно в следующем фрейме.
Оценка параметров модели осуществляется методом максимального правдоподобия. В R нет метода оценки GARCH-M модели, поэтому воспользуемся скриптом (кодом на R) профессора R.Tsay Чикагского университета, взятого с сайта https://faculty.chicagobooth.edu/ruey.tsay/teaching/introTS/. Оценим параметры GARCH-M процесса на примере рейтов RUB/USD. 10 - минутные данные января-февраля 2016 года
library(fGarch)
source('/ShurD/ALection/garchM.R')
data <- read.csv('c:/ShurD/Alection/RubleRates.csv',sep=';')
rubusd <- data$Rub.Usd
returns <- rubusd * 100
matplot(returns,type="l",lwd=3,col = "blue",main = "RUB/USD Rates")
g <- garchM(returns)## Maximized log-likehood: -248.9429
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## mu 0.11735896 0.09305440 1.26119 0.207242
## gamma -0.37890877 0.32345781 -1.17143 0.241425
## omega 0.00856316 0.00424741 2.01609 0.043790 *
## alpha 0.03897431 0.01574786 2.47489 0.013328 *
## beta 0.93515465 0.02128800 43.92873 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Оценка квадрата волатильности рубля
matplot(g$sigma.t,lwd=3,col = "blue",type="l",main ="Estimated process sigma_t^2")
Остатки после оценки модели
matplot(g$residuals ,lwd=3,col = "blue",type="l",main ="Estimated process sigma_t^2")
QQ-график остатков
qqnorm(g$residuals ,lwd=3,col = "blue")
Экспоненциальная GARCH модель
Другой моделью, позволяющей учитывывать гетероскедастичность в прогнозах это экспоненциальная GARCH модель, предложена Nelson (1991). Пусть \[u_{t}=\epsilon_{t}\sqrt{h_{t}}\]
где \(\epsilon_{t}\) -независимые одинаково распределенные случайные величины с нулевым средним и единичной дисперсией. Nelson предложил следующую модель для логарифма \(h_{t}\) \[ln h_{t}=\zeta+\sum_{j=1}^{\infty}\pi_{j}[|\epsilon_{t-j}|-E|\epsilon_{t-j}|+Λ\epsilon_{t-j}]\] здесь все \(\pi_{j}>0\). Параметр \(Λ\) играет следующую роль. Если \(Λ=0\) ,то вклад в общую волатильность проложительных значений \(e_{t-j}>0\) будет таким же, как и отрицательных \(e_{t-j}<0\). При \(-1 < Λ<0\) вклад положительных увеличивает волатильность меньше, чем от отрицательных. При \(Λ<-1\) приводит к тому, что положительные уменьшают волатильность, отрицательные увеличивают.
Естественной параметризацией данной модели является представление бесконечного полинома \(\pi(B)\) как отношение двух полиномов конечного порядка \[\pi(B)=((\alpha(L))/(1-\delta(L)))=((\alpha_{1}B+...+\alpha_{m}B^{m})/(1-\delta_{1}B-...-\delta_{r}B^{r})\]
аналогично тому, как это делалось в GARCH модели В этом случае исходную модель можно переписать следующим образом
\[ln h_{t}=k+\delta_{1}ln h_{t-1}+...+\delta_{r}ln h_{t-r}+ \alpha_{1}[|e_{t-1}|-E|e_{t-1}|+Λe_{t-1}]+ ...+\alpha_{m}[|e_{t-m}|-E|e_{t-m}|+Λe_{t-m}]\]
где \(κ=[1-\delta_1-...-\delta_{r}]\zeta\).
Параметры экспонециальной GARCH модели можно оценить методом максимального правдоподобия, естественно вид функции правдоподобия зависит от предположений относительно распределения шума \(e_{t}\). Нельсон, в качестве распределения для \(e_{t}\) предложил использовать обобщенное распределение ошибок, нормальзованное так, чтобы математическое ожидание было равно нулю, а дисперсия единице. \[f(x)=((νe^{-(1/2)|(x/λ)|^{ν}})/(2^{((ν+1)/ν)}λΓ((1/ν))))\]
здесь \(λ\)- константа равная
\[λ={((2^{(-2/ν)}Γ((1/ν)))/(Γ((3/ν))))}^{(1/2)}\]
При \(ν=2\) получим, что \(λ=1\). В этом случае \(f(x)\) - примет вид плотности стандартного нормального распределения. При \(ν<2\) плотность \(f(x)\) имеет более “толстый хвост” по сравнению со стандартным нормальным, а \(ν>2\) наоборот, более “тонкий” Для стандартного нормального распределения (случай \(ν=2\)) математическое ожидание для модуля ошибки имеет вид
\[E|e_{t}|=√((2/π))\]
а при произвольном \(ν>0\)
\[E|e_{t}|=((λ2^{(1/ν)}Γ((2/ν)))/(Γ((1/ν))))\]
Плотность обощенного распределения ошибок при параметре \(\nu=1,2,3\) выглядит следующим образом
x = seq(-4, 4, length = 501)
y = dged(x, mean = 0, sd = 1, nu = 1)
plot(x, y, type = "l", , col = "blue", main = "GED", ylab = "Density",
xlab = "")
y = dged(x, mean = 0, sd = 1, nu = 2)
lines(x, y, col = "red")
y = dged(x, mean = 0, sd = 1, nu = 3)
lines(x, y, col = "green")
abline(h = 0, col = "grey")
Для иллюстрации рассмотрим следующую модель Пусть моделью дневного дохода \(r_{t}\) по некоторой акции является следующая модель \[r_{t}=a+br_{t-1}+\delta h_{t}+u_{t}\] где ошибка \(u_{t}=e_{t}√(h_{t})\) и модель для \(h_{t}\) \[ln h_{t}-\zeta_{t} = \delta_1(ln h_{t-1}-\zeta)+\delta_2(ln h_{t-2}-\zeta)+ α_1(|e_{t-1}|-E|e_{t-1}|+Λe_{t-1})+ α_2(|e_{t-2}|-E|e_{t-2}|+Λe_{t-2})\]
\(\zeta\) и \(ρ\)- параметры, которые можно оценить методом максимального правдоподобия. Логарифм функции правдоподобия задается выражением \[L=T{ln((ν/λ))-(1+(1/ν))ln2-ln(Γ((1/ν)))}- (1/2)∑_{t=1}^{T}|(((r_{t}-a-br_{t-1}-\delta h_{t}))/(λ√(h_{t})))|^{ν}-(1/2)∑_{t=1}^{T}ln h_{t}\]
Где последовательность \(h_{t}, t=1,...,T\) получают итерационно
Для численного примера на R рассмотрим известный нам 10 минутные доходности курса рубля относительно доллара США в январе-феврале 2016 года
library(rugarch)
source('/ShurD/ALection/garchM.R')
data <- read.csv('c:/ShurD/Alection/RubleRates.csv',sep=';')
rubusd <- data$Rub.Usd
returns <- rubusd * 100
matplot(returns,type="l",lwd=3,col = "blue",main = "RUB/USD Rates")
Гистограмма рейтов по сравнению с плотностью нормального распределения
histPlot(as.timeSeries(returns))
говорит нам, что при использовании нормального распределения мы вряд ли получим лучшую модель в этом ниже и убедимся.
Оценим EGARCH(1,1) модель с AR(1) моделью для рейта. Для начала будем предполагать нормальность \(\epsilon_t\)
spec <- ugarchspec(mean.model = list(armaOrder = c(1, 0)),variance.model = list(model ="eGARCH",garchOrder = c(1,1)),distribution.model = "norm")
fit <- ugarchfit(returns, spec = spec)
fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : eGARCH(1,1)
## Mean Model : ARFIMA(1,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.052200 0.020117 2.5949 0.009463
## ar1 -0.092695 0.054395 -1.7041 0.088361
## omega 0.009938 0.005634 1.7642 0.077704
## alpha1 0.114235 0.021212 5.3854 0.000000
## beta1 1.000000 0.000145 6879.6670 0.000000
## gamma1 0.042805 0.010376 4.1254 0.000037
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.052200 0.024685 2.1147 0.034457
## ar1 -0.092695 0.045577 -2.0338 0.041971
## omega 0.009938 0.006941 1.4319 0.152186
## alpha1 0.114235 0.031558 3.6198 0.000295
## beta1 1.000000 0.000386 2590.3076 0.000000
## gamma1 0.042805 0.013376 3.2001 0.001374
##
## LogLikelihood : -243.1919
##
## Information Criteria
## ------------------------------------
##
## Akaike 1.5722
## Bayes 1.6433
## Shibata 1.5715
## Hannan-Quinn 1.6006
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1989 0.6556
## Lag[2*(p+q)+(p+q)-1][2] 0.6398 0.9217
## Lag[4*(p+q)+(p+q)-1][5] 1.3701 0.8771
## d.o.f=1
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.003803 0.9508
## Lag[2*(p+q)+(p+q)-1][5] 1.153276 0.8241
## Lag[4*(p+q)+(p+q)-1][9] 1.956446 0.9101
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 1.388 0.500 2.000 0.2387
## ARCH Lag[5] 1.683 1.440 1.667 0.5456
## ARCH Lag[7] 1.944 2.315 1.543 0.7292
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.6749
## Individual Statistics:
## mu 0.12625
## ar1 0.16489
## omega 0.11873
## alpha1 0.12247
## beta1 0.09609
## gamma1 0.12532
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.1993 0.2313
## Negative Sign Bias 0.0929 0.9260
## Positive Sign Bias 0.1230 0.9022
## Joint Effect 3.3602 0.3393
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 40.67 0.002673
## 2 30 50.67 0.007649
## 3 40 61.68 0.011802
## 4 50 79.69 0.003643
##
##
## Elapsed time : 0.481446Повторим оценку, но теперь будем использовать обобщенное распределение ошибок для \(\epsilon_t\)
spec1 <- ugarchspec(mean.model = list(armaOrder = c(1,0)),
variance.model = list(model ="eGARCH",
garchOrder = c(1,1)),
distribution.model = "ged")
fit1 <- ugarchfit(returns, spec = spec1)
fit1##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : eGARCH(1,1)
## Mean Model : ARFIMA(1,0,0)
## Distribution : ged
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.008340 0.002761 3.02116 0.002518
## ar1 -0.000509 0.002979 -0.17084 0.864352
## omega -0.002298 0.016336 -0.14070 0.888109
## alpha1 0.080526 0.018853 4.27127 0.000019
## beta1 0.999998 0.013865 72.12224 0.000000
## gamma1 0.083768 0.013395 6.25362 0.000000
## shape 1.018154 0.101977 9.98415 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.008340 0.000379 22.03140 0.000000
## ar1 -0.000509 0.000228 -2.22958 0.025775
## omega -0.002298 0.010975 -0.20942 0.834118
## alpha1 0.080526 0.020550 3.91857 0.000089
## beta1 0.999998 0.010145 98.57209 0.000000
## gamma1 0.083768 0.036043 2.32415 0.020117
## shape 1.018154 0.114797 8.86920 0.000000
##
## LogLikelihood : -221.8654
##
## Information Criteria
## ------------------------------------
##
## Akaike 1.4439
## Bayes 1.5270
## Shibata 1.4430
## Hannan-Quinn 1.4771
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.4342 0.5099
## Lag[2*(p+q)+(p+q)-1][2] 0.7848 0.8589
## Lag[4*(p+q)+(p+q)-1][5] 1.3386 0.8841
## d.o.f=1
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.01127 0.9155
## Lag[2*(p+q)+(p+q)-1][5] 0.92018 0.8775
## Lag[4*(p+q)+(p+q)-1][9] 1.47373 0.9577
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.9041 0.500 2.000 0.3417
## ARCH Lag[5] 1.2348 1.440 1.667 0.6648
## ARCH Lag[7] 1.3361 2.315 1.543 0.8536
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.9065
## Individual Statistics:
## mu 0.12582
## ar1 0.29797
## omega 0.09129
## alpha1 0.09235
## beta1 0.09233
## gamma1 0.19266
## shape 0.57881
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.9387 0.3486
## Negative Sign Bias 0.1150 0.9085
## Positive Sign Bias 0.4794 0.6320
## Joint Effect 2.6541 0.4481
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 22.75 0.2486
## 2 30 22.65 0.7922
## 3 40 35.68 0.6220
## 4 50 37.10 0.8938
##
##
## Elapsed time : 1.352539По критерию Акаике вторая модель предпочтительней первой модели(значение критерия меньше), хотя во второй модели больше параметров (добавился параметр shape \(\mu\)). Логарифм функции правдоподобия больше, что также в пользу второй модели.